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I introduce a modification of continuous matrix product states (CMPS) that makes them adapted to 
relativistic quantum field theories (QFT). These relativistic CMPS can be used to solve genuine 1+1 
dimensional QFT without UV cutoff and directly in the thermodynamic limit. The main idea is to 
work directly in the basis that diagonalizes the free part of the model considered, which allows to fit 
its short distance behavior exactly. This makes computations slightly less trivial than with standard 
CMPS. However, they remain feasible and I present all the steps needed for the optimization. The 
asymptotic cost as a function of the bond dimension remains the same as for standard CMPS. I 
illustrate the method on the self-interacting scalar field, a.k.a. the ¢3 model. Aside from providing 
unequaled precision in the continuum, the numerical results obtained are truly variational, and thus 


provide rigorous energy upper bounds. 


I. INTRODUCTION 


Tensor network states (TNS) have become the tool of 
choice for many-body quantum systems on the lattice 
[1, 2]. They provide a sparsely parameterized manifold 
of wavefunctions that is well suited for the approxima- 
tion of the low energy states of local Hamiltonians [3, 4]. 
The bond dimension D controls the size of the manifold 
and thus the expressiveness of the class of states. In 1+1 
spacetime dimensions, matrix product states (MPS) [5] 
are the simplest instance of TNS. They underpin the ear- 
lier density matrix renormalization group [6-8] and pro- 
vide one of the most efficient numerical methods to study 
quantum chains. Their extension for 2 + 1 dimensional 
problems, the projected entangled pair states (PEPS) [9], 
while more difficult to optimize, have been used success- 
fully in a wide range of non-trivial instances including 
the Hubbard model. Aside from their numerical use, ten- 
sor network states are powerful theoretical tools [10], e.g. 
for the classification of topological |11, 12] and symmetry 
protected [13-16] phases of matter. 

In light of the progress they enabled in the discrete, 
it is tempting to use TNS to solve problems in the con- 
tinuum, and attack quantum field theories (QFTs). This 
can be done in two ways, (i) by discretizing the QFT first 
to then use the standard lattice TNS toolbox and finally 
extrapolate the results (ii) by taking the continuum limit 
of TNS first to apply them to the continuum model with- 
out extrapolation needed. The first option has already 
provided numerically accurate results in many quantum 
field theories [17-22]. However, working directly in the 
continuum offers crucial advantages: spatial symmetries 
are preserved and the perilous continuum extrapolations 
of the results are avoided. It is this neater second option, 
which consists in defining a manifold of states directly in 
the continuum, that I am interested in here. 

In 2010, Verstraete and Cirac introduced continuous 
matrix product states (CMPSs) [23], which are a proper 
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continuum limit of MPS particularly adapted to 1 + 1 
dimensional non-relativistic field theories (like the Lieb- 
Liniger model). The corresponding numerical toolbox 
was then developed and extended in the past ten years 
[24-26]. More recently, an extension to d+ 1 (d > 2) di- 
mensions pushing PEPS to the continuum, was put for- 
ward in [27], but without a general numerical toolbox 
yet. 

The previous continuum ansatz are adapted to non- 
relativistic theories only, and cannot be used for rela- 
tivistic QFTs without an additional momentum cutoff 
A. There are many ways to understand this limitation 
which I summarize in sec. II B. For relativistic QFTs, this 
partially defeats the purpose of going to the continuum 
in the first place, since one still needs to extrapolate the 
final results as A —> +00. This is particularly disappoint- 
ing given that, at least in 1 + 1 dimensions, relativistic 
QFTs are straightforward to define “all the way down”, 
without any cutoff [28, 29]. 

My main objective in this paper is to introduce a mod- 
ification of CMPS, the relativistic CMPS (RCMPS), that 
is adapted to relativistic QFT in 1 + 1 dimension, and 
does not require any additional cutoff, UV or IR. In a 
nutshell it is a state |Q, R} parameterized by two D x D 
complex matrices Q, R and defined as 


IQ, R) = tr {Pexp | [gert roat] io). 


I will explain this formula in more detail later. Let 
me just mention here that a'(x) is the Fourier trans- 
form of the momentum creation operator al diagonal- 
izing the free part of the relativistic Hamiltonian un- 
der consideration and |0), the associated Fock vacuum. 
Crucially, a(x) is not the creation operator 7 (x), local 
in the canonically conjugated fields, and used for non- 
relativistic CMPS. Using a instead of w is the only dif- 
ference between CMPS and RCMPS. We will see later 
why using a is a good idea and the consequences this 
choice has. 

I can mention one advantage already. Because 
RCMPSs allow to work directly with the true theory 


without cutoff or extrapolation, the results one obtains 
are genuinely variational, and thus provide rigorous en- 
ergy upper bounds. Further, while results based on a 
discretization or with a finite UV cutoff have a limited 
validity at high momenta, the RCMPS are exact in the 
UV. 


The price to pay is that RCMPS are obtained from 
CMPS via change of basis w, Wt — a,a' that is not local. 
In the new basis, the Hamiltonian density of a relativistic 
QFT is no longer strictly local, and only exponentially 
decaying. In practice, the lack of strict locality of the 
Hamiltonian density introduces minor complications in 
the computations. However, and perhaps surprisingly, 
these complications do not change the asymptotic scaling 
of the cost compared to regular CMPS, which remains 
x D. Further, the energy density seems to converge 
fast as a function of D as one would expect [30], while 
the cost which means RCMPS provide an efficient class 
of states for relativistic QFTs. 


In addition to their direct use for relativistic QFT in 
1+ 1 dimensions, demonstrated in this paper, RCMPS 
could be used as an auxiliary step to solve non-relativistic 
QFTs in 2+ 1 dimensions. Indeed, continuous general- 
izations of PEPS currently lack a general purpose opti- 
mization algorithm. The missing routine is a function to 
solve a relativistic QFT in one dimension less, that is in 
1+ 1 space-time dimensions [27]. The relativistic nature 
of this auxiliary theory comes from the need for Euclidean 
invariance in the 2 space dimensions of the original non- 
relativistic theory [27]. For that task, RCMPS would fit 
the bill better than CMPS and preserve exact Euclidean 
invariance at short distances. 


The present paper is structured as follow. I first dis- 
cuss standard CMPS and their limitations for relativistic 
theories in sec. II before introducing and discussing the 
RCMPS ansatz in sec. III. I first explain how to com- 
pute expectation values in sec. IV and then show how 
to optimize the state to find ground states in sec. V. 
The numerical results for the ground energy and cor- 
relation functions, which result from this optimization, 
are presented for the ¢* model in sec. VI. I finally dis- 
cuss some extensions in sec. VII, namely a slight varia- 
tion of RCMPS with more general basis, and explain how 
one could obtain more general observables by extending 
known CMPS techniques. I conclude with a general dis- 
cussion of the method in sec. VIII, comparing it with 
renormalized Hamiltonian truncation (RHT). A quicker 
presentation of the results, emphasizing the importance 
and difficulty of the variational method, is presented in 
a companion letter [31]. 


II. STANDARD CONTINUOUS MATRIX 
PRODUCT STATES 


A. Definition and main properties 


Continuous matrix product states (CMPS) were intro- 
duced by Verstraete and Cirac [23] in 2010. They extend 
the successful matrix product states ansatz from the lat- 
tice to the continuum. Concretely, for a bosonic quan- 
tum field theory on the closed line [—L, L] with (non- 
relativistic) creation and annihilation operators (t, y), 
a CMPS is the state 

0)» 


(2) 
where by definition [y(x),~'(y)] = (x — y), P is the 
path ordering operator, and |0),, is the Fock vacuum an- 
nihilated by all the y(x). The state is parameterized by 
two D x D matrices Q and R, where D is the bond di- 
mension, and which contain all the variational freedom. 
The trace is taken over the associated D dimensional 
Hilbert space, and implements periodic boundary condi- 
tions (more general boundary conditions are possible by 
inserting a matrix B in the trace). 

CMPS can be derived from a genuine continuum limit 
of MPS and thus share many of their properties [32]. In 
particular, all the correlation functions of local operators 
can be computed explicitly and efficiently. This is seen 
by introducing the generating functional Z; j: 


L 
IQ, R) nr = tr {Pes I drQQ1+R8Qyİ(z) 
=i 


(Q, Ri exp (f 7’ Yt) exp (S jY) |Q, R)nr 
(Q, RIQ, BR) ar i 


which can be used to compute all normal-ordered corre- 
lation N-point functions of local operators, e.g. 


(Q, Ridt(x)b(y)|Q,R) ne 
(Q, RIQ, R) ur 
ô ô 
KOLOK 


2515 = (3) 


(ilayh) = 


(4) 


J,j/=0 


Using Wick’s theorem, one can show [32] that this gen- 
erating functional has an exact expression: 


L 
Zj j=tr {Pel | dzT+ j(xz)R@1+ j'(x)1 8 Rl 
=L 
(5) 


where T = Q@14+1@Q+R® Ris the transfer operator 
and the trace is now taken over the tensor product of two 
copies of the original D dimensional Hilbert space. For 
example, for the simple two-point function of eq. (4), 
using the exact expression gives, for x > y: 


(yi (x)b(y)) =tr [e-ta @ Rje ITR Q Ljet +T : 
(6) 


The state |Q, R)nr is parameterized in a redundant way, 
in the sense that different Q, R pairs give the same state. 
In particular, conjugating both matrices with an invert- 
ible matrix U keeps the state invariant. This freedom 
can be used to choose a gauge in which only the anti- 
Hermitian part K of Q is free [33]: 


Q := -iK - SRIR. (7) 


In this gauge, the transfer operator T is of the Lindblad 
form, hence negative, and generically has a single largest 
eigenvalue Ag = 0 with associated left and right eigenvec- 
tors (€9| and |ro). Using this form provides convenient ex- 
pressions in the thermodynamic limit, as e”™ > |r9)(0o| 
when L — +00. For example, eq. (6) simplifies to: 


(wt (x) W(y)) = Woll S Re? (R®DIro). (8) 


In what follows, and in particular for the upcoming rela- 
tivistic extension, I will always work directly in the ther- 
modynamic limit. 

A CMPS is typically used to find the ground state of 
a non-relativistic QFT. The archetypal example is the 
Lieb-Liniger model with Hamiltonian Hu: 


a [ Abad t eytip- uit (9) 


= f tn. (10) 


With the generating functional, it is possible to compute 
the expectation value of the 3 terms of the Hamiltonian 
density (hLL)Q,r = f(Q, R), where f is a simple function 
containing expectation values of products of Q and R as 
in eq. (8). One then minimizes the energy density over 
the matrices Q (or K) and R to get close to the ground 
state: 


|ground) ~ |Q, R) where Q, R = argmin (hLL)Q,R- 

(11) 
The approximation can get arbitrarily good as the bond 
dimension (and thus the size of the variational manifold) 
is increased because CMPS are dense in the Fock space 
[32]. Further since the method is variational, we also 
know that the energy we find for finite D always upper 
bounds the true ground energy 


E0 := (hit) ground < (han) o,R- (12) 


In practice, one can simply feed the expression of the 
energy density as a function of Q and R to a standard 
numerical minimizer as was done originally in [23]. How- 
ever it is typically much more efficient, especially for large 
D, to use more elaborate tangent space methods [34]. I 
explain how they work in sec. V. In a nutshell, the lat- 
ter require computing the exact gradient with respect to 
the variational parameters as well as the natural metric 
on the tangent space of CMPS induced by the real part 
of the Hilbert space scalar product. The optimization is 


then done through gradient descent on the correspond- 
ing differentiable manifold. At this stage, all that matters 
for us is that it can be done and that it is efficient: one 
can easily optimize CMPSs of large bond dimensions [24] 
(without barren plateaus or other pathology). 


B. UV problems in relativistic theories 


CMPS are well suited for non-relativistic theories, which 
they approximate well “all the way down”, without short 
distance cut-off. However, once we move to relativistic 
QFT, it becomes necessary to introduce a UV cut-off. 
This is best understood on the free boson, which was 
discussed in the CMPS context in [35]. The free boson 
Hamiltonian is 


Hm = F T? + (p) + me, (13) 
2 JR 


where m is the mass and 7, ¢ are canonically conjugated 
(d(x), r(y)] = iô(x— y). To deal with such a Hamiltonian 
with a CMPS, it is tempting to express the field opera- 
tor and its conjugate as a function of a non-relativistic 
creation-annihilation pair yt, w: 


b= awry), (14) 
r= $w- (15) 


which introduces a new arbitrary mass scale A. In this 
basis, the Hamiltonian Hp is still the integral of a local 
density, and thus it is straightforward to evaluate the 
energy density on a CMPS |Q, R) nr. 

This unfortunately leads to mild and serious diver- 
gences. The first mild one is that Hp, is a priori not nor- 
mal ordered in the ~',w operator basis. This is simply 
solved by considering : Hmw:y instead. The serious diver- 
gence comes from the fact that : Hg: not only contains 
the standard non-relativistic kinetic energy 0,~'0,~ but 
also 0,WO0,W + h.c.. The latter is divergent when eval- 
uated on a generic CMPS [33]. This second divergence 
can be cured by adding an adapted UV regulator to the 
Hamiltonian and considering [35] 


1 
Hi, = sf” + (Org)? +m? do? + — (Oer), (16) 
regulator 


which kills the problematic terms. To reach a fixed pre- 
cision, one then needs to increase the bond dimension 
D as the cutoff A is lifted. Similarly for fermionic theo- 
ries, it was observed in [36] that one needs to add a UV 
cutoff to the Hamiltonian, this time not to ensure the 
finiteness of the results, but to make the optimization 
problem well behaved. In a nutshell, without a cut-off, 
the CMPS reduces the energy density by approximating 
larger and larger momentum modes as the optimization 


proceeds, completely missing the IR which contributes 
to the energy only in a subleading way. Independently of 
CMPS, this sensitivity to high frequencies in variational 
methods was considered by Feynman as one of the main 
difficulties preventing its use in relativistic QFT [37]. 

Zooming out from the CMPS peculiarities, this sit- 
uation is not surprising, since we are working in the 
wrong operator basis. Indeed, the operator : Hm:y is not 
bounded from below with this specific ordering. Hence 
even if we manage to have a finite energy density when 
computing CMPS expectation values, we are always in- 
finitely far from the true ground state. Further, the free 
boson is a conformal field theory (CFT) at short dis- 
tances (the massless free boson). It is thus no surprise 
that a CMPS, which is adapted to systems with a gap and 
exponentially decaying correlation functions, completely 
fails to capture the short distance behavior of relativistic 
models. Yet another way to see the problem is that the 
ground state of the free boson has an infinite density of 
non-relativistic particles, whereas a CMPS always gives 
a finite density. 

Note that these UV problems are not made easier if one 
adds a relevant interaction, as the latter becomes negligi- 
ble at short distances. The relativistic continuous matrix 
product state ansatz will solve these short distance prob- 
lems already present in the free theory, but also allow to 
deal with interactions without additional difficulty. 


III. RELATIVISTIC CONTINUOUS MATRIX 
PRODUCT STATES 


A. Intuition 


In the standard textbook approach to quantum field the- 
ory [38], the problem of infinite particle density is solved 
at the very beginning by changing of basis (and in fact 
even of Hilbert space). Typically, one expands the field 
operator in new creation-annihilation modes that diago- 
nalize the Hamiltonian: 


= 1 J 1 ika —ikx Tt 

d(x) = On fa on (c an +e al) (17) 
= 1 Wk ikx _ piks (ft 

n(£) = z [yt (e ape al.) , (18) 


where wp = Vk? +m? and [ax,al,] = 27rô(k — k’). In 
condensed matter physics, this would be called a Bogoli- 
ubov transform. In this new basis, the Hamiltonian is 
diagonal 


1 i i 
Ha = ae f drer akak + Akap (19) 
oT R 2 


Then, normal ordering Hp, —: Hfp:a removes the infi- 
nite vacuum energy contribution and one finds that the 
ground state is simply the Fock vacuum annihilated by 
all the a,’s: |ground) = |0),. The Hamiltonian : Hfp:a 


further is a well defined operator on the Fock space gen- 
erated by creating relativistic particles from the vacuum. 

The crucial observation is that in 1+1 dimensions, this 
normal ordering procedure is typically! sufficient to cure 
all the divergences that can appear, even when adding 
interactions. In particular, the ¢* Hamiltonian H which 
we will consider in more detail later, 


H =: Hoiatg | pta, (20) 
R 


is a perfectly legitimate and regular Hamiltonian on the 
free Fock space. From now on, I will omit the subscript 
on the normal ordering which will always be done with 
respect to a,a unless otherwise stated. 

The Fock space generated by acting with al on |0)q is 
much more adapted to relativistic theories than the Fock 
space generated by acting with yt on |0),, because the 
former solves the scale invariant short distance behavior 
of the theory exactly with its vacuum. However, the 
operators az, al are not adapted to CMPS because the 
latter encode the state in real space, and not momentum 
space. Further, the Hamiltonian H is not translation 
invariant in momentum, which is the situation CMPS are 
most adapted for. A tempting workaround is to simply 
Fourier transform a, to get a(x) 


1 ik 
a(x) = xfa etap, (21) 


which verifies [a(x),at(y)] = 6(a — y). Note the cru- 
cial fact that y(x) 4 a(x), as the Fourier transform (21) 
does not contain the factors wg. Intuitively, at(x)|0}a 
corresponds to a (bare) relativistic particle localized in 
x. Working with this new operator basis is the key to 
make CMPS adapted to relativistic theories. 


B. Definition and basic properties 


The previous discussion naturally leads us to introduce 
the relativistic CMPS (RCMPS) ansatz: 


IQ, R) = tr f Pexp [fagen +Real(o)} 0)a 
(22) 


where at (x) is the Fourier transform of the creation oper- 
ator diagonalizing the free part of the relativistic Hamil- 
tonian under consideration. Note that we choose to work 
directly in the thermodynamic limit, and thus the state 
IQ, R} has no cut-off, UV or IR (although the latter is 
trivial to reintroduce). 

The properties of the state are the same as before if 
one replaces ~ by a, as these operators obey the same 


1 Normal ordering is sufficient to make polynomials in the field 
well defined. For more general potentials, like cos(b¢), normal 
ordering is sufficient only for b small enough. 


algebra. A large part of the theory of CMPS can thus be 
reused. In particular, all normal-ordered N-point corre- 
lation functions of a, a can be computed efficiently using 
the same generating functional. Additionally, we have 
that RCMPS are dense in the appropriate QFT Fock 
space (the manifold is maximally expressive) and thus 
that the precision can be arbitrarily refined by increas- 
ing D. 


C. Consequence on the Hamiltonian density 


The main difference with the non-relativistic case is that, 
once written as a function of a(x), the Hamiltonian of 
a (massive) relativistic theory is not strictly local, but 
only exponentially decaying. Indeed a relativistic Hamil- 
tonian is local in the field ¢(x) and its conjugate, which 
do not have a local expression as a function of a(x). More 
precisely, 


= 1 1 ika —ikx + 
da) = = fayz (c ak +e al) 


_ = ae (c—Maly) + e#H@-Wal (y)) 
= f ay Hey) faly) + ah(y) (23) 


where J(z) is a smooth kernel for x 4 0 (and not a Dirac 
distribution) which I will make explicit later. 

Consequently, the Hamiltonian density of a relativistic 
theory with terms of order up to n in the field @ can 
be written as a function of a(x) with n integrals, a bit 
informally: 


H= fa h(x) (24) 


= fè avr denK (er ak tn)a) (¢1)...a\P (an) 
(25) 


where al) is a compact notation to convey the fact that 
both creation and annihilation occur in a sum, and K is 
a kernel that decays exponentially as a function of the 
difference of its arguments. 

Is this non-locality a problem? Technically, it certainly 
induces complications in evaluating the energy density: 
we have exact expressions for N point functions of a, at 
which we then have to integrate over some kernel (instead 
of taking them at equal point for a local density). This 
is however not insurmountable, and as I will show in 
the next section, one can still compute the Hamiltonian 
density at a cost x D3. 

More crucially, and provided we can optimize them, do 
we expect the expressive power of RCMPS to be lower for 
such Hamiltonians? There is no reason to think so, and in 
fact CMPS have already be applied to Hamiltonians with 
exponentially decaying interactions in the non-relativistic 
context [39]. Intuitively, we are merely introducing a new 


length scale m~! that replaces the lattice scale in lattice 
models. This is more physical than introducing a much 
smaller and arbitrary cut-off scale A~' < m7! as was 
done previously. Other choices of length-scales that could 


further improve precision are discussed in sec. VII. 


IV. EFFICIENTLY COMPUTING RCMPS 
EXPECTATION VALUES 


A. Naive direct evaluation 


The operators a'(x),a(y) have the same commutation 
relations as the (x), Y(y) of standard non-relativistic 
CMPS and thus all the formulas follow. In particular, 
one can straightforwardly compute the exact expectation 
values of normal-ordered products of at and a using the 
generating functional (3). 

Obtaining simple functions of the field ¢, e.g. expec- 
tation values of normal-ordered monomials (:¢”:) is less 
straightforward because the field ¢ is obtained from a, a! 
with a convolution (23). Hence, to compute the expecta- 
tion value of (:6”:), one a priori needs to compute n inte- 
grals (in fact n — 1 using translation invariance) of exact 
aat correlation functions (a“)(2)a\ (x2) ---a (ap)). 
This is feasible for low bond dimension and can be used 
as a sanity check, but is prohibitively expensive for a full 
optimization of the state’. 


B. Vertex operators 


For efficient computations of functions of the field ¢, the 
starting point is to compute expectation values of vertex 
operators: 


(Vs) := (Q, R| :e°?): |Q, R), (26) 


which we may evaluate in x = 0 without loss of gener- 
ality because of translation invariance. Such an expec- 
tation values seems even more difficult to compute than 
field monomials at first sight. Indeed, using the naive 
approach above and expanding the exponential (26) pro- 
vides n integrals at each order n. However, it turns out 
one can directly compute the expectation value without 
expanding the exponential. 

To this end, we simply express the vertex operators as 
a function of a(x) 


POR = / T / ee Mala) + e'ttat(a)|; 
— b | dda (| ae b J di H(a)a(x) 
(27 


) 


2 This was, unfortunately, the strategy I first followed. 


where J is a real function 


1 dk _ i, 
= — An 2 
Jle) 2T QW. (28) 
Ky /4(|2/m)) 


= DAAT EA fap P 
and K,(x) (not to be confused with the matrix K) is the 
modified Bessel function of the second kind. We see from 
(27) the crucial fact that the expectation value of a vertex 
operators is simply the generating functional itself, with 
bJ as source, i.e. (V,) = ZyzoJ, thus explicitly 


(Vo) = tr [Pep f aT +0s(0 (Rer+1eÑ)] ; 


(30) 
Inside the trace, this is simply the solution of an ordinary 
differential equation that can be solved numerically. 

To reduce the computational cost, one can use the 
standard trick of matrix product states which consists 
in mapping the tensor product Hilbert space CP @ CP 
to the space of matrices p acting on CP. Introducing the 
super-operators 


1 
L- p= —i[K, p] + RpR' — 5 (R'Rp+pR'R) (31) 


R(x) - p = J(x)(Rp + pR') (32) 


we have 
(Va) =tr{Pexp | [ave + ora] “push (33) 


where pss is the stationary state of £ normalized with 
trace 1. Rewriting the path-ordered exponential as the 
solution of an ODE we have 


(Vi) = lim tr [pa] (34) 


where limy-+—oo Px = pss and 


gale 5 E Pa + OR(2) > pa (35) 
The limit (34) is well defined because J(x), and thus 
R(x), decrease exponentially fast at infinity and are in- 
tegrable in 0. Because of this fast decay at infinity, one 
could in fact use any density matrix as initial state. Us- 
ing a simple ODE solver, e.g. a backward differential 
formula (BDF) solver, one can obtain the limit in (34) to 
arbitrary precision with only a reasonable number of sub- 
divisions. The total computational cost is proportional 
to the cost of applying the super-operator L + bR(x) on 
a density matrix and thus scales x D? only. 


C. Field monomials 


To compute field monomials, one can differentiate vertex 
operators with respect to their exponent b 


n 


a 
CO") = g (36) 


This allows to obtain (: ¢”:) by directly differentiating 
the ODE (35). Doing so yields 


(:¢":) = lim tr [e] (37) 


x£r— +00 


where p‘*) := OF p|,-9 obey n coupled matrix ODE 


d 3 
TI =L EARE) pe? (38) 
x 

with the convention that po) = pss and for k > 0, 


pE = 0. Solving the ODE above numerically provides 


arbitrarily accurate approximations of (: 6” :} at a cost 
xn x D?. 


D. Kinetic term 


In addition to exponentials and polynomials of the field 
ġ, it is important to be able to compute the expectation 
value of the free part of the Hamiltonian. 

For convenience, we consider directly the Hamiltonian 
for the massive free boson, but since the mass term 
m?(: ¢? :) is also computable, it could be subsequently 
subtracted to obtain the pure kinetic term. The free 
boson Hamiltonian can be expressed as a function of 
the momentum space creation and annihilation operators 
al, ak and reads 


1 
: Ap := — faros alap. (39) 
2T 


The corresponding Hamiltonian density hp,(x) is 


 Agy(a) := ~ faray wp e®Y-Mat(ya(x). (40) 


As before, we would like to write this density as a deriva- 
tive of a vertex operator. If we try to mirror the reasoning 
of the previous subsection, we face the issue that the nat- 
ural source J that appears now is the Fourier transform 
of y/wp. This is not a function but only a distribution. 
To get a true function, we can divide and multiply by 
w? = m? + k?°, and interpret the k? term on the numer- 
ator as a 0,0, derivative. This gives 


1 dk ria 
ihola) = 5 f dy E 0-9 (m? + 0,8,)al (yale). 


(41) 
We are back to an expression that depends on the source 
J(x) that we introduced before (28) 


(hy 3) =2m? ( f des(e)al(a)f autojat) Y) 


42 ( ‘| E E J ays(u)2ya(a)) = 


except this time derivatives of a,a' appear as well. This 
gives 


(43) 


where J; ; is the generating functional of normal ordered 
correlation functions of zat (x), 0ya(y). This generating 
functional also has an exact expression, which is easily 
derived by differentiating correlation functions obtained 
from Z with respect to position 


Y; j=tr {Pexo( [7 +Q, R] @1+j7'1® 2.) \ 

| (44) 
As before, Be Bey Zoi J,baJ and 3 on J,boJ Can be ob- 
tained by solving simple ODEs. To this end, we intro- 
duce px := p?!°2 with initial condition p_. = pss and 
dynamics 


d 
— pr = L - pe tb I(x) Rpg + b2J (x£) pz R (45) 


da 
and Oy := gbibe with initial solution o_. = pss and 
dynamics 
d 


de? = L-Os +b J(x)[Q, R]os +b2J(x)o2[R', Q'] (46) 


We further introduce notations for the partial deriva- 
tives p10) = Ob; Plbs,2=0; pO) = bə P|b1,2=0 and 
p™® := 8b, Ab. Plo, 2—0. They obey 


d 

p09) = L. pl + Ja) Roo (47) 
T 

d 

a —f. plod site J(x)poR' (48) 
d 

4 LD = L. pD + sayy” + slap RE (49) 


dx 


The same system of ODEs can be obtained for a, replac- 
ing R by [Q, R] and Rt by [R', Qt]. Finally, the expecta- 
tion value we are looking for is obtained from the trace 
of the solutions 

(: hp :) =2 lim tr [mog +D]. (50) 

xz— +00 

Hence the expectation value of the massive free boson 
Hamiltonian density (: hg, :) can be computed by solving 
2 systems of 3 coupled matrix ODEs, and thus can be 
obtained to arbitrary precision at a cost x D®. 


V. OPTIMIZATION 
A. Failure of naive optimization 


Using the results in the previous section, one obtains 
an expression for the energy density of the form (h) = 


f(Q, R) where f is a function of the matrices R and Q (in 
practice K) that can be evaluated efficiently on a classi- 
cal computer at a cost x D?. One may thus simply input 
this function to a standard minimizer, that will typically 
use a gradient computed by finite differences and hope 
for the best. This is what was done in the original paper 
on CMPS applied to the Lieb-Liniger model [23]. For our 
model, this approach works reasonably well up to D = 4, 
after which all the standard optimizers (e.g. L-BFGS or 
conjugate gradient) get stuck in plateaus. To understand 
why it happens and go beyond such small values of D, 
we need to do better, and use a tangent space approach. 


B. Tangent space approach 


It is well known that the notion of steepest descent in 
optimization depends on a choice of metric. More pre- 
cisely, if we want to minimize a function f(x) where 
z = {a }"_, is a vector of parameters (think of the co- 
efficients of R,Q), we can go down the steepest descent 
direction u defined by 


u = argmin (V f, U) - (51) 


lul=1 


This scalar product on the tangent space (u,v), = 
Guv(z)u"v” and associated metric g,, are a priori ar- 
bitrary. The notion of “steep” depends on a metric, and 
what is steep for the “right” metric g,, may look like a 
plateau for the naive 4,,, metric if gv is singular. If there 
are many parameters, the naive metric has no reason to 
be good. 

But what is the right metric? It is one where the 
distance between parameter values is proportional to 
how much they change the function one optimizes. An 
excellent choice of metric is thus given by the Hes- 
sian Hess,, := 0,0, f of the function one is optimiz- 
ing. Taking this metric gives the descent direction u œ 
—[Hess~']#”0, f where [Hess~']#”Hess,) = ôb, which 
corresponds to the famous Newton method. This ma- 
trix is costly to estimate for RCMPS, because it requires 
computing « D* derivatives of the energy density, in- 
stead of x D? if we only compute the gradient. 

There is another natural option that comes from the 
fact that, in our case, the tangent space is also a Hilbert 
space [40]. Indeed, let us write |x) = |Q, R} a state in 
the manifold of RCMPS. Then a natural tangent space 
metric is simply the Hilbert one 


Gy (T) := Re [(A,(2|)(A.|2))] . (52) 


It provides a notion of distance between parameter val- 
ues associated to how much they change the quantum 
state (instead of the energy). Further, in our case, it 
can be computed straightforwardly (it is instantaneous 
in comparison with the computation of the gradient). 

A more physical justification for the use of this metric 
is that it corresponds to (approximate) imaginary time 


evolution [40], which converges exponentially fast for a 
gapped system and an expressive enough state manifold. 
Indeed, upon an infinitesimal imaginary time evolution 
dr an RCMPS |z} evolves into |”) —drH|x). This latter 
state no longer belongs to the RCMPS manifold, and to 
get an approximate evolution we need to project down 
the evolution to the tangent space. More precisely, we 
want to find a direction u € R in the tangent space 
such that u0,,|z) ~ —H |z}. It is obtained by minimizing 
u”ð le} + H|zx)||?, which gives 

ul = —|g7]"” ð ( (x| H |x)) (53) 
provided g,,, is invertible which we will assume here. 
This projected imaginary time evolution corresponds 
to the (imaginary) time dependent variational principle 
(TDVP) in the tensor network context [34]. In my opin- 
ion, the advantage of seeing imaginary TDPV simply as 
gradient descent with a different metric is that it makes 
it obvious the time step does not need to be infinitesimal, 
and can be chosen optimally with a line search. 

In practice, I observed for D < 4 that quasi-Newton 
methods, which try to estimate the best metric (the Hes- 
sian) from the gradient at different iterations, are still 
reasonably efficient. For larger D, I found that the metric 
Juv becomes very singular near the ground state, which 
may explain why quasi-Newton methods fail to estimate 
the Hessian (which is likely very singular as well) and 
get stuck in plateaus. However, as we will see in V E, the 
tangent space approach I presented converges fast even 
for large values of D as one would expect. For the opti- 
mization RCMPS in moderately large D, it is thus better 
to have an exact “good” metric, than an approximation 
of the best metric. 


C. Computing the metric 


The metric can be computed easily following [34]. The 
first step is to define the tangent space vectors 


IQ, R). 
(54) 


ô 
lV, Whar = [ax i + Wa 


ô 
5Qaa(2) f 5Rag(a) 


The complex matrices V,W parameterize the direction 
in the tangent space, and Q, R the point on the RCMPS 
manifold. We work in the translation invariant case, 
where Q, R are taken position independent at the end, 
but the position argument x in (54) is necessary to know 
how operators are ordered. A crucial fact is that the tan- 
gent space is overparameterized and, with the left canon- 
ical choice (7) we took for Q, ie. Q = —iK — RİR, 
one is free to fix V = —R'W without losing a linearly 
independent direction [34]. We may thus drop V as a 
parameter, as it is fixed by W. With this choice, one can 
show [34] that the overlap between tangent vectors takes 


the particularly simple form 


(Wi|W2)a,r = (€0|W2 ® Wi|ro) 
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= tr[W2pssW]] (a 
where pss is the (normalized) stationary state of the Lind- 
bladian £ defined in (31). We thus have 2D? directions 
on the tangent space corresponding to the real and imag- 
inary parts of the coefficients Wag. The metric is simply 
the bilinear map taking two W and outputing the real 
part of (55). Note that the metric depends on the state 
only and is cheap to compute. Indeed, it does not re- 
quire the resolution of ODEs which are the numerical 
bottleneck of RCMPS. 


D. Computing the gradient with an adjoint method 


To compute the gradient of the energy density in the 2D? 
independent directions, one a priori needs « 2D? com- 
putations of expectation values each with a cost x D3. 
However, using standard adjoint methods (a.k.a. back- 
propagation), one can compute the complete gradient 
with the same asymptotic cost as computing the energy, 
hence x D3. In principle, this could be done by using 
a complex ODE solver compatible with automatic differ- 
entiation. In practice, because the ODEs involved have 
a special form, it is easy, efficient, and illuminating to 
implement the adjoint method directly. I illustrate the 
idea on the example of the computation of the gradient 
of a vertex operator, but it applies equally easily to the 
gradients of field monomials and kinetic term. 

The expectation value of a vertex operator on a 
RCMPS is 


(Vi)a.r = tr {Pexp ji dz £ + PRC) “puch. (56) 


Let us consider the gradient in the W direction 
Vw(Vo)Q,r Which is defined implicitly via 


(Vi)o+ev.ntew = (Vodo, R +eVw(Vb)Q,r + Ole). (57) 


Differentiating directly (56) yields 


+00 nb Yy b 
Vw (V) = /dytr {Pel Tayt Pel ps} 


(58) 
with the notation L’(x) = L + bR(x) and 


1 
Vwl?(y)-p=Vo+pVvi+ z (RoW! + WpR') (59) 


+ bJ(y) (Wp + pW") , 


recalling again that V = —R'W. This gradient of £°(y) 
appears in (58) between two evolution super-operators. 
Because a trace is taken at the end, we can replace the 


last part of the evolution from y to +00 by its adjoint 
acting on the identity 


Vel = forf [peh a] 


<vw£%ty): [Pela pa] (60) 


where the adjoint £°*(y) of £°(y) is defined as 


1 
L™(y)-O = Qi0+0Q+ ROR + bJ(y) [RIO + OR]. 
(61) 
a b 
Writing as before py = Pel bee Pss the solution of the 


+00 nbs 
forward problem and O, = Bel. ae 1 the solution of 
the backward problem we get 


Vi) = fayte[O,Vw0"u)- ry] (62) 


This can be further simplified exploiting the expression 
of VwL?(y) (59), and one obtains all the components of 
the gradient from 2 matrices My and My: 


Vw (Vo) = tr [MwW + My: W'] (63) 


where 


1 
Mw =| dy — pyOyR' + shy RO, + bI(y) pyOy 
i (64) 
My =f dy —ROypy + zO PP + bI(y)Oy py - 


In practice, one computes py and Oy by solving the 
corresponding ODEs. The matrices My and My are 
then obtained by evaluating the integrals in (64) with 
an efficient numerical method like the tanh-sinh quadra- 
ture [41]. This gives an algorithm with a cost x D? to 
compute the full gradient of the expectation value of a 
vertex operator. The gradients of the kinetic term and 
of other potentials can be computed efficiently with the 
same method. 


E. Algorithm 


We now have all the pieces to understand the optimiza- 
tion algorithm. The first step is to start from an initial 
guess. One can certainly do much smarter, but I started 
from uniformly random R and K matrices. This gives 
a very high starting energy density, but it fortunately 
decreases fast enough that initialization is a secondary 
concern for the bond dimensions I probed. 

The second step is to compute the descent direction, 
obtained by acting with the inverse metric on the gradi- 
ent. The gradient is computed with the backpropagation 
method described before, and is the most costly step, 
while the computation of the inverse metric is an imme- 
diate algebraic operation using (55). 


The third step is to move R and Q in the descent di- 
rection. The step size need not be small, as in imaginary 
time evolution, and it is chosen to approximately yield 
the maximal energy decrease at each step. In practice, 
I used a backtracking line search with Armijo-Goldstein 
condition to find the the right amount to move at each 
step. Note that this is very similar to what was done 
by Ganahl et al. in [24] for standard CMPS and the 
Lieb-Liniger model. Like them, I observed that this ap- 
proach speeds up the optimization by roughly 2 orders of 
magnitude compared to imaginary time evolution with 
an optimal but fixed time step. Typically, results are 
converged after ~ 10? — 104 iterations depending on the 
coupling (convergence is slower near criticality), bond di- 
mension (convergence is slower for large D), and random 
initial seed. 


VI. APPLICATION TO THE 
SELF-INTERACTING SCALAR 


A. The model 


To assess the soundness of RCMPSs, we consider the 
simplest non-trivial QFT in 1 + 1 dimensions, the self- 
interacting scalar field with Hamiltonian 


_ [fm , O) m 
n= | |+ 5 +H +99" F (65) 


The normal ordering is again done with respect to the 
creation and annihilation operators al, ap which diag- 
onalize the quadratic part of the Hamiltonian and are 
defined in (17). 

This model is a good case study because it is simple 
to define, even rigorously, as H is a genuine renormalized 
Hamiltonian (self-adjoint, finite energy density). Yet, the 
model is not integrable, and carrying accurate computa- 
tions out of the perturbative regime is non-trivial. The 
self-interacting scalar has been studied with a wide va- 
riety of methods: renormalized Hamiltonian truncation 
(without space-time discretization but finite size) [42], 
infinite matrix product states (with space discretization 
but no finite size cutoff) [17, 43], Monte-Carlo (with 
space-time discretization and finite size) [44-46], ten- 
sor network renormalization (with space-time discretiza- 
tion and finite size) [21, 22], and, of course, (resummed) 
perturbative expansions (without cutoff, but perturba- 
tive) [47]. 

Out of these works, let us mention two that are par- 
ticularly relevant for our study as they are carried in the 
Hamiltonian formalism. The study of Milsted et al. [17] 
is the the closest, in terms of method used, to what we 
shall do: the authors discretize the model in space, and 
find the ground state with translation invariant matrix 
product states, thus without IR cutoff, reaching unbeat- 
able precision at the time. The drawback is the need to 
extrapolate the continuum limit, i.e. the UV cutoff. On 


the other hand, the remarkably pedagogical Hamiltonian 
truncation study of Rychkov and Vitale [42] is carried di- 
rectly in the continuum with the Hamiltonian (65). In 
a nutshell, the authors introduce an IR and an energy 
cutoff that make the Hilbert space finite, and exactly 
diagonalize the resulting finite Hamiltonian matrix. In 
addition, they introduce a smart renormalization proce- 
dure which makes the convergence of the results faster 
as the cutoffs are lifted. The drawback is that the size 
of the Hilbert space (and thus the cost) is exponential 
in both cutoffs (I discuss it further in VIII). Further, 
while the renormalization procedure and careful extrap- 
olations drastically improve precision, they destroy the 
variational nature of the results. 

My objective is not so much to improve upon these 
studies in terms of raw numerical precision than in terms 
of conceptual simplicity, robustness, and scaling: with 
RCMPS one can in principle find the ground state of 
(65) directly in the continuum limit, without UV or IR 
cutoff, and in a variational way (that is, with rigorous 
energy upper bounds). 


B. Ground energy density 


The ground state energy density is finite and negative 
for the model we consider. This is clear given that the 
Fock vacuum for the free part |0),, which is no longer the 
ground state for g # 0, already gives a zero energy be- 
cause the Hamiltonian is normal ordered. Optimizing the 
RCMPS for m = 1 indeed shows that the energy density 
is negative and decreases as the coupling is increased (see 
Fig. 1), at first quadratically in the coupling constant g 
(as expected from perturbation theory [42]). 

Importantly, the convergence in bond dimension is fast 
for all values of the coupling, even deep in the non- 
perturbative regime (which kicks in roughly at g > 0.1). 
As we can see in Fig. 1 the energy as a function of g is 
already qualitatively correct for D = 5, and the points 
at larger bond dimensions (D = 10, 15, 20) are essentially 
indistinguishable on this plot. At large coupling (g > 3 
the results are substantially below those obtained with 
renormalized Hamiltonian truncation (RHT) [42], which 
means they are more accurate since the method is truly 
variational. 

To estimate the error, it is not possible to compare with 
an exact solution (#3 is not integrable), nor with earlier 
numerical results, e.g. RHT, which are less precise even 
in their latest high precision development [48]. To get 
an accurate point of comparison, I simply considered a 
large D estimate of the energy density as reference to es- 
timate the error at lower bond dimensions. For D = 32, 
I obtained the rigorous bounds (h) =, < —0.039354 and 
(h) g=2 < —0.157214, which can be used as fairly good es- 
timates of the true values. The resulting errors for lower 
bond dimensions are shown in Fig. 2 for g = 1 and g = 2 
and provide hints that the convergence to the ground en- 
ergy is close to exponential in the bond dimension, or 
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FIG. 1. Ground state energy density of the f model as a 
function of the coupling g for m = 1. The RCMPS results 
are compared with the renormalized Hamiltonian truncation 
calculations from [42], obtained with an IR cutoff L = 10. 
Points are sampled more densely around the critical coupling 
Qc & 2.77. 


at least faster than any power law as one expects from 
the discrete [30]. Retrospectively, this fast convergence 
implies that our D = 32 points of comparison are ex- 
act enough, with an error to the true ground state much 
smaller than that of the lower D points whose error is 
estimated. 

Note that the energy density obtained with RCMPS is 
the renormalized one and thus would be very difficult to 
estimate with a similar precision with lattice methods. 
Indeed, the latter give access to the total energy density, 
which diverges in the continuum limit. One would need 
to subtract the diverging tadpole part to get the finite 
renormalized contribution. As one gets closer to the con- 
tinuum limit, obtaining this finite correction to a fixed 
precision requires a prohibitively high relative precision 
on the total energy. 

Finally, note that a second order phase transition oc- 
curs for ge ~ 2.77 (for m = 1) [22]. Just like Rychkov and 
Vitale [42] with RHT, and as expected from the similar- 
ity with the 2d Ising model, we see no sign of this tran- 
sition in the ground state energy density. As we will see, 
the transition appears more clearly once we look at ob- 
servables like the expectation value of the field ¢ (which 
behaves like the Ising magnetization). 


C. Observables 


Once an RCMPS is optimized to approximate the ground 
state of a given Hamiltonian, expectation values of oper- 
ators come essentially for free. As an illustration, I show 
the expectation values of ¢ (Fig. 3) and : ¢?: (Fig. 4) 
but one could equally easily consider : 647: or :cosh(¢):. 

The expectation value of ¢, which is the equivalent of 


relative error for (h) 


5 10 15 20 25 
D 


FIG. 2. Relative error in the energy density as a function 
of the bond dimension D. The error is computed by using 
as reference the results at D = 32: (h)g=1 < —0.0393547 
and (h)g=2 < —0.157214. The latter should have a substan- 
tially lower error (~ 107° and œ 1074 respectively) and thus 
be exact enough for the approximate computation of relative 
errors. 
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FIG. 3. Absolute value of (¢) taken in the approximate 
RCMPS ground state. The symmetry breaking as the cou- 
pling is increased is manifest. Away from the critical point, 
convergence is extremely fast as a function of the bond di- 
mension. 


the Ising magnetization for the ¢+ model, is the most in- 
structive. It shows a clear spontaneous symmetry break- 
ing around the expected coupling ge ~ 2.77. To locate 
it more precisely (with a precision closer to the lattice 
extrapolations [22]) and estimate the critical exponents, 
one would need to compute more data points near the 
critical point, for a wide range of D, and use finite en- 
tanglement scaling techniques [35, 43, 49]. This is left 
for future work. Consequently, the exceptional numeri- 
cal accuracy of the ansatz is so far limited to the gapped 
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FIG. 4. (:¢°:) in the approximate RCMPS ground state. The 
phase transition can also be seen, although less clearly than 
in Fig. 3, through the divergence of the first derivative near 
criticality. 


phases on both sides of the critical point. 


VII. EXTENSIONS 
A. Adjustable characteristic length-scale 


The core idea of RCMPS compared to CMPS is to 
use creation operators such that the theory is exactly 
solved at short distance. We used the pair a(x), at (x) 
associated to the free part of the theory with mass m. 
As we discussed, this introduces a length-scale m~1 for 
the exponentially decaying support of the Hamiltonian 
density, which is reminiscent of the lattice scale for lattice 
models. But this scale is arbitrary in our case and we 
could have chosen a different pair of creation-annihilation 
operators as long as the UV behavior is still exact. 

The simplest extension is to consider mass-adjustable 
RCMPS, that is states of the form 


Q Rm) =tr{Pexp | [igen ro ahto] Yoa 


(66) 
and where am are the operators diagonalizing the free 
theory of mass ñ, and |0}m is the associated Fock vac- 
uum. More explicitly, the field operators can be ex- 
panded in this new operator basis as 


_ 1 1 ike _ ika tî 

olx) = =j DA (e am. k +e a a (67) 
— 1 Wh ikx | ika Tt 

n(x) = J > (c am,k e az k) ; (68) 


where w(k) = vk? + M?. At short distances, large k, 
the state |Q, R, Mm) still solves the theory exactly just like 


the RCMPS |Q, R) = |Q, R, m), but the variable mass or 
inverse length-scale m 4 m gives an additional degree of 
freedom one can optimize to better fit the IR. 

The computations with the mass variable RCMPS are 
more difficult than with the standard RCMPS, and the 
complications mainly come from the fact that the Hamil- 
tonian density is not normal-ordered for the operators 
am, a}. Evaluating the Hamiltonian density is thus te- 
dious but doable, and I could carry the optimization. 
However, for the few values of coupling I tried, I obtained 
rather underwhelming results, only marginally improving 
the precision at the cost of a substantial increase in com- 
plexity and slower optimization. 

A more thorough study should be done in future works. 
One could expect that carefully optimizing the length- 
scale ñm! near criticality would yield more meaningful 
improvements, since the gap becomes much smaller than 
the mass m appearing in the Hamiltonian (corresponding 
to the one-loop renormalized mass). A more systematic 
exploration of Bogoliubov transformations would be in- 
teresting as well. In principle, one can change the a, at 
into new operators b,b! by replacing wg in the mode ex- 
pansion (17) with any function 0, with the same large 
k behavior so that expectation values are still UV finite. 
This should provide a substantial gain in expressiveness 
at fixed bond dimension, but the computations would be 
more involved. 

Finally, I would like to emphasize that this idea of non- 
local change of basis to make the continuum limit well 
behaved or even simply increase expressiveness, which 
comes at the cost of making the Hamiltonian exponen- 
tially decaying instead of local, could be used on the lat- 
tice as well. 


B. Excitation spectrum and beyond 


Once optimized, a RCMPS can be used to compute all 
correlation functions at equal time for no additional min- 
imization cost. As I argued, this also gives some dynam- 
ical information because of Lorentz invariance, but can 
one get more? 

In principle, one can use TDVP in real-time to evolve 
states which means we have access to all dynamical prop- 
erties. Note however in that case that one cannot use the 
trick of taking optimally large time steps, since, in real 
time, we no longer have a global minimization problem. 

Using standard tangent-space CMPS techniques, one 
also has access to the excitation spectrum [34]. In a nut- 
shell, the idea is to diagonalize the Hamiltonian in the 
tangent space |V,W)g,r which is a vector space orthog- 
onal to the ground state (in the gauge we chose). This 
would give approximations only to the excited states with 
zero momentum, but simple local modifications allow to 
target states of non-zero momenta as well [34]. Carrying 
such computations requires evaluating matrix elements of 
the form (V',W"\h|V,W) gr, which should be doable. A 
comparison of this spectral data from the one one could 
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extract from the two-point function would be interesting. 


C. Other quantum field theories 


I illustrated the use of RCMPS on the self-interacting 
scalar field only, but it can in principle be applied to 
almost all theories in 1 + 1 dimensions. Other bosonic 
theories with polynomial interactions (e.g. : 6° :) can 
be treated without new technique. Scalar theories with 
exponential potentials, like the Sine-Gordon and Sinh- 
Gordon models can also be dealt with immediately since 
expectation values of vertex operators are straightfor- 
ward to compute. 

Fermionic theories could also be dealt with directly, 
with a minor subtlety related to regularity conditions 
[32], that is conditions on R, Q one has to impose to make 
expectation values finite. The Gross-Neveu model [50], 
which has already been studied with various UV cutoffs 
with tensor networks [36, 51], would be an interesting 
candidate to probe the behavior of a “just renormaliz- 
able” theory. Alternatively, a large class of Fermionic 
models, like the Thirring model, can already be dealt 
using bosonization and the results of the present paper. 

There are however serious difficulties remaining to ex- 
tend the method to relativistic QFT in 2+ 1 and 3+1 
dimensions. The first is specific to the Hamiltonian for- 
malism, as renormalized Hamiltonians are more difficult 
to define in higher dimensions: normal ordering is no 
longer sufficient and the renormalized Hamiltonian no 
longer acts on the free Fock space [52]. Nonetheless, 
recent progress was made using renormalized Hamilto- 
nian truncation [53] and lightcone conformal truncation 
[54], with promising numerical results, showing that the 
Hamiltonian approach is still reasonable in 2+ 1 dimen- 
sions and in the continuum. 

The second difficulty is related to continuous tensor 
network states themselves, that is the higher dimension 
equivalent of CMPS. In the non-relativistic setting, these 
states have been proposed in [27], but evaluating cor- 
relation functions is so far efficient only for Gaussian 
states [55]. Indeed, the equivalent of the finite trans- 
fer matrix T of CMPS in 2+1 dimensions is an operator 
acting on (two copies of) the Fock space of a 1+1 dimen- 
sional relativistic QFT. In principle, one could bootstrap 
the present approach and evaluate correlation functions 
(or the energy density) in 2+1 using a boundary RCMPS 
approach. In a nutshell, one would find the stationary 
state of the transfer matrix as a (large bond dimension) 
RCMPS. Every evaluation of an expectation value in 2+1 
would be done at the cost of a full optimization in 1+ 1. 
This is the typical dimensional reduction obtained with 
tensor network approaches, where a physical dimension 
is traded for a variational optimization. Clearly, optimiz- 
ing RCMPS is so far orders of magnitudes too slow to be 
used as a routine to be called many times in a full opti- 
mization process. I hope the present work can stimulate 
work in this direction. 


VIII. DISCUSSION 


In this paper, I have presented a new class of states, 
the relativistic continuous matrix product states, that is 
adapted to relativistic quantum field theories. It is usable 
directly in the thermodynamic limit (no IR cutoff) and is 
exact at short distances, all the way down (no UV cutoff). 
The bond dimension D controls the expressiveness of this 
class, and a state has 2D? independent parameters. 

The comparison between RCMPS and Hamiltonian 
truncation is a good illustration of the difference between 
linear and non-linear methods. Hamiltonian truncation 
(up to its renormalization refinements) is a linear ap- 
proach. The candidate ground state is expanded in a 
truncated Hilbert space. The energy is quadratic in the 
coefficients, and thus minimized efficiently as a linear 
problem. The price to pay is a lack of extensivity, the 
size of the Hilbert space grows exponentially as the size 
of the system increases for a fixed energy cutoff (and thus 
fixed precision). A side effect is that the number of pa- 
rameters needs to be huge to reach good precision, with 
the candidate ground state written as a linear superpo- 
sition of typically 10* to 10° states [48]. In contrast, 
using a continuous tensor network ansatz as we did, we 
work with a manifold of states, and the energy density 
is a highly non-linear function to optimize. Nonetheless, 
it can still be done efficiently as I showed in V. Further, 
having a non-linear class of states allows us to have an ex- 
tensive ansatz, where going to the thermodynamic limit 
does not increase the number of parameters (in the trans- 
lation invariant case). The results at D = 5 provide an 
approximation to the ground state with only 2D? = 50 
independent real parameters which is already more accu- 
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